function [OutDS, OutOptionGrid] = OutsideOptionModelSimple(DS,opt)
% creates outside option return and partial derivatives wrt price using the
% principle: outreturn = (local corn share)xCorn_yield + (local soy
% share)xSoy_yield.
% classification of fields is based on two variables: (local corn
% share)xCorn_yield and (local soy share)xSoy_yield


% Load data

%load([opt.DataFolder opt.MicroData])
load([opt.DataFolder opt.MacroData])

% Load long price series:
prices = [];
load([opt.DataFolder opt.PricesData])

% Variables from Agregate data to keep:
VarAgregKeep = {'munic', 'Year', 'micro_reg', 'meso_reg', 'uf', ... % Identifiers
                        'p_sugarbr_lag', 'p_sugarbr', 'sb1', 'sb12', 'sb1brl', 'sb12brl', ... % Sugar price data
                        'sb1_lag', 'sb12_lag', 'sb1brl_lag', 'sb12brl_lag', ... 
                        'c1brl_lag', 'c1brl', 'p_cornbr', ... % Corn prices
                        's1brl_lag',  's1brl', 'p_soybr', ... % Soy prices
                        'area_planted', 'corn_planted', 'soy_planted', 'area_farms_cs06', ...
                        'cs_corn_qty', 'cs_corn_area', ...  % Census 2006 corn info
                        'cs_soy_qty', 'cs_soy_area'}; % Census 2006 soy info}; % Cane, corn, soy and total farm land
                        
Agregate.Properties.VarNames{3} = 'Year';
Agregate.Properties.VarNames{1} = 'munic';

Agregate = Agregate(:,VarAgregKeep);

Agregate = sortrows(Agregate,{'munic','Year'});
Agregate = Agregate(Agregate.Year == 2006,:); % Use only census year (should not matter)

% Assign regional variables to DS
RegionCorrespondence = Agregate(:,{'munic','micro_reg','meso_reg','uf'});
[~,Ridx] = unique(RegionCorrespondence(:,1));
RegionCorrespondence = RegionCorrespondence(Ridx,:);

DS = join(DS,RegionCorrespondence,'LeftKeys','munic','RightKeys','munic');

% Regions for outside option return:
RegionVariable = opt.RegionVariable;
RegionsAll = double(Agregate(:,RegionVariable));
RegionsUni = unique(RegionsAll);
nr = length(RegionsUni);

Agregate.AgregArea = nan(size(Agregate,1),2);

for r = RegionsUni'
    Agregate.AgregArea(RegionsAll == r,1) = nansum(Agregate.cs_corn_area(RegionsAll == r));
    Agregate.AgregArea(RegionsAll == r,2) = nansum(Agregate.cs_soy_area(RegionsAll == r));
end

% compute relative shares:
Agregate.AgregShare = Agregate.AgregArea./(sum(Agregate.AgregArea,2)*ones(1,2));
Agregate.AgregShare(sum(Agregate.AgregArea,2) == 0,:) = zeros(sum(sum(Agregate.AgregArea,2) == 0),2);

% n = size(Agregate.AgregShare,1);
% Agregate.AgregShare = [zeros(n,1),ones(n,1)];

DS = join(DS,Agregate(:,{'munic' 'AgregShare'}),'LeftKeys','munic','RightKeys','munic');

DS.SY(:,1) = DS.AgregShare(:,1).*DS.gaez_aecol_mze/1000;
DS.SY(:,2) = DS.AgregShare(:,2).*DS.gaez_aecol_soy/1000;

% Classify fields into basereturn (separately for corn and soy):
NumPartialStates = opt.NumPartialStates;
LimClass = quantile(DS.SY,linspace(0,1,NumPartialStates+1)');
[PartialIdx,SYMidPoint] = ClassStates(DS.SY,LimClass);

Ns = size(SYMidPoint,2);
Ny = NumPartialStates;
PartialStatesMap = MultIdx(Ny,Ns);

OutIdx = ones(size(DS,1),1);
for s = 1:Ns
    OutIdx = OutIdx + (PartialIdx(:,end+1-s)-1).*( Ny^( Ns - s ) );
end

OutOptionGrid.Partialcorn = SYMidPoint(PartialStatesMap(:,1),1);
OutOptionGrid.Partialsoy  = SYMidPoint(PartialStatesMap(:,2),2);

AgriPrices = double(prices(:,{opt.CornPrice, opt.SoyPrice}));

pavcorn = nanmean(AgriPrices(:,1));
pavsoy  = nanmean(AgriPrices(:,2));

OutOptionGrid.Base = pavcorn*OutOptionGrid.Partialcorn + pavsoy*OutOptionGrid.Partialsoy;
OutOptionGrid.PriceAverage(1) = pavcorn;
OutOptionGrid.PriceAverage(2) = pavsoy;

OutDS.Index = OutIdx;
OutDS.BaseReturn = DS.SY(:,1)*pavcorn + DS.SY(:,2)*pavsoy;

function MIdx = MultIdx(Ns,Ny)
    Nstar = Ns^Ny;
    MIdx = NaN(Nstar,Ny);
    UVndx = (1:Ns)';
    MIdx(:,1) = repmat(UVndx, Ns^(Ny-1), 1);

    for j = 2 : Ny
       n1          = Ns ^ (j - 1);
       n2          = Ns ^ (Ny - j);
       MIdx(:,j)  = repmat(kron(UVndx, ones(n1,1)), n2, 1);
    end
end

end

